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Adaptive  Wavelet  Galerkin  Methods 
on  Distorted  Domains: 

Setup  of  the  Algebraic  System 


Stefano  Berrone  and  Karsten  Urban 


Abstract.  We  use  the  algorithm  of  Bertoluzza,  Canuto  and  Urban  [2] 
for  computing  integrals  of  products  (of  derivatives)  of  wavelets  in  order  to 
solve  elliptic  PDEs  on  2D  distorted  domains.  We  construct  a variant  of  the 
original  method  which  turns  out  to  be  more  efficient.  Several  numerical 
results  are  presented. 


§1.  Introduction 

Adaptive  wavelet  Galerkin  schemes  have  quite  recently  been  proven  to  offer 
great  potential  for  numerically  solving  boundary  value  problems  for  partial 
differential  equations.  On  the  one  hand,  strong  analytical  properties  such 
as  convergence  and  optimal  efficiency  have  been  proven  for  elliptic  operators 
[6,9].  On  the  other  hand,  first  numerical  tests  also  on  non-tensor  product 
domains  indicate  the  applicability  of  such  methods,  [1], 

However,  the  major  obstacle  so  far  is  the  efficient  computation  of  the 
entries  of  the  stiffness  matrix  and  the  right-hand  side  of  the  corresponding 
algebraic  systems.  In  fact,  it  turns  out  that  these  entries  are  more  expensive  to 
compute  than,  e.g.,  in  the  case  of  adaptive  Finite  Element  Methods.  In  [2],  a 
method  to  adaptively  approximate  and  compute  these  entries  was  introduced 
and  analyzed;  numerical  results  were  given  for  a ID  example.  In  this  paper, 
we  study  the  application  of  the  algorithm  in  [2]  for  2D  ‘distorted’  domains, 
which  are  parametric  images  of  the  unit  square.  This  allows  the  study  of  the 
influence  of  ‘realistic’  parametrizations  of  non-tensor  product  domains  on  the 
assembling  of  the  algebraic  system.  We  incorporate  some  improvements  over 
the  original  method  in  [2]  to  increase  efficiency,  and  present  various  numerical 
results. 
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§2.  Adaptive  Approximation  of  the  Algebraic  System 

Given  a linear  boundary  value  problem  in  a bounded  Lipschitz  domain  Cl  C 
Hn  (n  > 1),  its  numerical  approximation  by  a variational  method  (Galerkin, 
Petrov-Galerkin,  weighted  residuals,  ...)  requires  the  computation  of  integrals 
of  the  form 


f bafiix)  Dau(x)  D^v{x)dx 
J n ’ 


where  a,  [3  6 IN’1  are  suitable  multi-indices,  D°  = 


/ fp{x)  D^v(x)  dx, 
Jn 

0IWI 


(1) 


with  llall  := 


Bxf1  ■ ■ ■ dx°" 

Old han,  u and  v are  suitable  trial  and  test  functions  belonging  to 

and  respectively,  ba tp  e L°°( fi)  and  fp  6 L2(Cl). 

We  consider  a wavelet  Galerkin  method  with  trial  and  test  spaces  S a = 
span  generated  by  adaptively  choosing  a finite  subset  'Pa  = {ip\  : A € A} 
within  a wavelet  basis  'P  = {ip\  : A € J } in  L2(Q),  i.e. , A C J (see,  e.g., 
[5,10]).  The  wavelets  are  assumed  to  have  the  appropriate  regularity  for  the 
above  integrals  to  be  well  defined. 

The  construction  of  such  wavelet  bases  on  fairly  general  domains  Q,  is  not  a 
trivial  task.  However,  quite  recently  significant  progress  has  been  made  on  this 
topic,  see  [3,4,8,12]  and  also  [13]  for  a somewhat  different  approach.  The  main 
idea  behind  all  the  constructions  in  the  first  cited  papers  is  domain  decom- 
position and  matching.  The  domain  f!  is  subdivided  into  N non-overlapping 
subdomains  f2j.  Each  subdomain  is  mapped  to  the  n-dimensional  reference 
cube  f2  :=  [0,  l]n  by  means  of  smooth  parametric  mappings 

Fj  : A — » Hi,  & = #$),  Gi-F'1.  (2) 

Then,  each  A € J,  restricted  to  fl,  is  the  image  through  F,  of  a linear 
combination  of  tensor  product  wavelets  Tpx  on  Cl,  i.e.,  if  A = {j,  k)  ( j =:  |A| 
denoting  the  level  and  k the  location  in  space  as  well  as  the  type  of  wavelet), 
then 

A'eS(i.A) 

where  S(i,  A)  is  a suitable  set  of  indices  of  the  form  A'  = (j,  k’)  with  k'  € Cl 
and  jy  . are  suitable  coefficients  independent  of  j (see  e.g.  [3,4]). 


2.1.  Reduction  to  univariate  integrals 

We  will  only  consider  the  calculation  of  the  integral  on  the  left-hand  side  of 
(1)  which  enters  into  the  stiffness  matrix.  The  entries  for  the  right-hand  side 
are  treated  analogously,  [2].  Hence,  replacing  u by  i})\  and  v by  Vv  on  the 
left-hand  side  in  (1)  for  some  A,  p € A,  we  get 

f N 

oa,m  :=  / batf3(x)Daipx(x)D0i)ll{x)dx  = ^ Tv,t  7£',« 

*=1  a'eT(i,a),  A'eS(t.A), 

/3'£T(i,/3)  ft'eS(i,ii ) 

X / ba>p(Fi(x))  dai  {x)  dpi  (x)  |JFj(i)|  Da  ipy{x)  D0  %j)p(x)dx,  (3) 

J n 
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where  the  sets  T(i,  a),  T(i,  (i)  are  defined  by  the  chain  rule,  and  da>  and  dp, 
are  smooth  functions  depending  on  Gi  and  its  derivatives.  The  integrals  which 
appear  on  the  right-hand  side  of  (3)  take  the  form 

/ c(x)Da$^(x)DP$p(x)dx,  where  ^(x)  = JJ »%.(£,)  (4) 

i=i 

and  are  univariate  scaling  functions  or  wavelets  on  [0,1].  Now,  we  use 
the  Two-Scale-Relation  for  the  wavelets  to  express  them  in  terms  of  scaling 
functions  on  the  next  higher  level,  i.e., 

^a(£)  = Y (5) 

A6AX 

where  m-x  ^ are  the  refinement  coefficients.  Here,  the  index  set  C Zi^|+1  is 
determined  by  the  Two-Scale-Relation  and  lj  denotes  the  set  of  all  scaling 
function  indices  on  a level  j.  Hence,  d~x  . becomes 

4,£  = Y Y mA,AmA,A  / c{x)b*<pk{x)bP<pp(x)<lx.  (6) 

AeAx  AeA^  JSi 

The  computation  of  each  integral  on  the  right-hand  side  of  (6)  would 
be  highly  efficient  if  we  could  reduce  it  to  a product  of  univariate  integrals, 
but,  in  general,  the  function  c is  not  a tensor  product  of  univariate  functions. 
However,  we  can  expand  it  in  an  appropriate  tensor  product  wavelet  basis 

0*  :=  {et  : §l(x)  = f[  (7) 

i= 1 

(where  jdt.  are  again  univariate  scaling  functions  and  wavelets  on  [0, 1],  re- 
spectively, possibly  different  from  ,0^)  as  follows 

e(£)  = Y c^o(£)-  (8) 

i>e  j- 

Then,  we  approximate  c locally  on  :=  supp  (p-x  fl  supp  (pp  by  a finite 

sum  Qa*  c,  obtained  by  restricting  the  sum  in  (8)  to  a finite  index  set  A*  C 3* 
(depending  on  c as  well  as  on  a,  f3,  A,  p,  A and  p),  whose  precise  definition 
will  be  given  below.  Correspondingly,  d~x  - is  approximated  by 

:=  'Y  Y i Q\-c{x)D&Vx{x)D0<pp{x)dx,  (9) 

A€AxA£Aa  Jn 
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which  is  a finite  linear  combination  of  products  of  univariate  integrals  of  the 
following  form 


J fa)  iO^ixi)  i ) ( &i ) dii , 


1)  • ■ 


(10) 


An  algorithm  for  computing  such  integrals  can  be  found  in  [2],  However,  here 
we  use  biorthogonal  B-spline  wavelets  [7,11]  as  trial  and  test  functions.  Due 
to  their  explicit  representation,  efficient  direct  formulas  for  the  integrals  (10) 
are  available  and  have  been  used  for  the  subsequent  numerical  experiments. 

Let  us  mention  that  the  above  strategy  slightly  differs  from  [2]  since 
here  we  approximate  c locally  on  Sx  whereas  in  [2]  this  is  done  on  the 

somewhat  larger  domain  ^ :=  supp  tpx  H supp  . This  new  method  ensures 

automatically  that  only  non-zero  integrals  are  computed,  avoiding  a wide 
number  of  checks,  which  explains  why  the  present  method  is  more  efficient 
than  the  original  one. 


2.2.  Adaptive  approximation  of  the  stiffness  matrix 

Now,  we  are  going  to  describe  the  construction  of  the  index  set  A*  introduced 
above.  To  this  end,  we  have  to  introduce  some  notation.  Let  us  set 

I{\,  A)  :=  {v  € J*  ■ |supp  §1  n supp  <pxr\  supp  | > 0} , (11) 

and  j :=  min{|A|,|A|}  as  well  as  J :=  max{|A|,  |£l}-  Let  Ro  be  the  number 
of  zero  moments  of  [5,10].  Moreover,  let  Tx  and  To  be  the  largest  integers 
such  that  (px  € iyTi’°°(S7)  and  6*0  € respectively.  Then,  we  set 

R:=  min{^,TA-||d||,T^-||4||}. 

We  make  the  following 

Assumption  1.  The  system  0*  defined  in  (7)  allows  the  characterization  of 
the  Besov  space  Bqq(fl)  for  indices  (a,  q)  in  a certain  range  Sq-  C M+  x (0, 1], 
i.e.,  the  Besov  seminorm  | ■ \B„  ^ has  the  representation 

/ \ 1/9 

l*k,(A)  ~ f E 2li>|n(9/2-1)  M* ) . * € BZJfl).  (12) 

The  following  notation  will  be  frequently  used  in  the  sequel.  For  l = 
1,  ...,£/,  we  consider  (possibly  different)  wavelet  bases  (ty  = {ftp \t : X(  € tJ). 
Then,  for  X(  6 iJ,  i — 1, . . . , L,  we  define 


i(A, At)  :=  I *'  ifinLiMPP*,l>0, 

[ 0,  otherwise. 
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Now,  we  are  in  a position  to  define  the  set  A*.  Let  us  fix,  once  and  for  all, 
independently  of  A and  ft,  a non-increasing  i1  (Wo)-sequence  S — (6e)eeiNo 
with  strictly  positive  elements,  whose  ^1(IVo)-norm  is  close  to  1.  Let  c € 
for  some  (cr,q)  € <S©*.  For  an  index  v S J* , we  define  its  relevance 

for  the  computation  of  Co  Op(x)  Daip^(x)  D^tp^x)  dx  as 


pff\ v)  :=j(A,A,»>)2",i>k,2-l£>ln^2-1)  \c>\'-<> 

(y-{R+n/2){\v\—J)  nnJ/2  olAKII^II-O  olAKIlAll-*)*-! 

■■  “ ~ iiAi-Iaii 

Finally,  for  any  e > 0,  we  define 

A*  :=  {u  € J(A,A)  = #Aa  #Am)  or  H < J}, 


(13) 


which  concludes  the  construction  of  an  adaptive  approximation  d7)^  of  d\4 L. 

Remark  2.  The  construction  of  A*  according  to  (13)  seems  to  require  the 
explicit  knowledge  of  all  the  (infinite)  coefEcients  c-c  of  c.  However,  one  can 

estimate  a priori  a level  Jc  such  that  < £/(mx,\  mA.A  #Aa  #aa)>  ^ 

\v\  > Je.  Following  [2]  it  is  easy  to  show  that  this  is  valid  for 


Je~ 




i?  + cr  + n(l-i) 


(14) 


where  we  have  set 

Qe  :=  |log2e/(m^mA^#A^#AA)|  + (R  + n)J  + |A|(|d|  - s) 

+ IA|(I/3|  - s)  + log2(|c|^(ft)d-i_|M||)  + log2  Const. 

Replacing  d\ltt  in  the  computation  of  a\tll  in  (3)  by  d*X  jl  results  in  an 
adaptive  approximation  a\  of  a^iM.  As  already  mentioned,  one  can  construct 
an  adaptive  approximation  f)  for  the  entry  of  the  right-hand  side  f\  := 
ffi  fpix)  D13^ x(x)  dx  in  the  same  way. 

2.3.  Error  estimates 

Let  us  assume  that  the  boundary  value  problem  we  aim  at  approximating 
is  elliptic  of  order  2s.  Let  H£(Q)  be  the  closed  subspace  of  Hs(fl)  which 
accounts  for  the  given  boundary  conditions.  The  wavelet  basis  ’3>  intro- 
duced above  is  assumed  to  form  a Riesz  basis  of  this  space.  The  wavelet 
Galerkin  approximation  of  our  problem  is  obtained  by  replacing  by 

Sa  ■==  span{V>A  : A € A),  where  again  A is  an  adaptively  chosen  subset  of  J . 
The  corresponding  Galerkin  solution  will  be  denoted  by  «a  :=  '“aV’a- 

The  vector  Ua  :=  (ma)acA  is  obtained  by  solving  the  linear  algebraic  system 
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PhylfcaJ  Demur*  ft  ScAifon 


Fig.  1.  Physical  domain  (left)  and  exact  solution  (right)  for  the  numerical  tests. 


A,\Ua  = fA,  which  is  defined  in  a straightforward  manner.  Let  the  integrals 
which  appear  in  the  stiffness  matrix  as  well  as  the  right-hand  side  be  com- 
puted in  an  approximate  way,  as  described  above.  Denote  the  resulting  matrix 
by  Aa  and  the  resulting  vector  by  fj(;  let  be  the  solution  of  the  modified 
linear  system  Aau)(  = f)(  and  let  u*A  — £)AgA  u\ 

The  following  estimate  on  the  effect  of  the  described  approximation  of 
the  stiffness  matrix  and  the  right-hand  side  has  been  established  in  [2].  It  is 
readily  seen  that  it  also  holds  for  our  variant  of  the  original  method  in  [2]. 

Theorem  3.  Under  the  above  and  similar  assumptions  for  computing  the 
right-hand  side,  there  exists  e0  > 0 such  that  for  all  0 < e < £q-' 


lMA  ^aIs.P  <- 
ImaU.o  ~ 


(15) 


§3.  Numerical  Results 

In  this  section,  we  present  our  numerical  results.  We  consider  the  Poisson 
problem  with  homogeneous  Dirichlet  boundary  conditions  on  a domain  ft 
which  is  the  parametric  image  of  ft  under  a suitable  transformation.  The 
domain  is  displayed  in  Figure  1,  left.  The  boundary  of  ft  consists  of  two 
straight  lines  and  two  curved  parts.  We  computed  the  parametrization  of  the 
four  parts  of  the  boundary  and  then  the  parametric  mapping  F : ft  — > ft  is 
determined  by  transfinite  interpolation,  [14]. 

We  constructed  a solution  u which  satisfies  the  boundary  conditions  and 
which  has  a strong  layer  near  the  upper  right  corner  of  the  domain.  This 
function  is  shown  on  the  right  in  Figure  1.  Since  we  have  an  explicit  formula 
for  u,  we  determined  the  right-hand  side  f by  using  MAPLE  V. 

The  choice  of  these  parameters  allows  us  to  test  an  interesting  situation. 
Indeed,  the  parametric  mapping  is  obviously  far  from  being  a tensor  product. 
Hence,  we  can  study  the  influence  of  a ‘realistic’  transformation.  Even  though 
this  influence  was  studied  in  [2]  in  ID,  we  face  here  a non-tensor  product 
situation  for  the  first  time.  Moreover,  for  computing  the  right-hand  side,  two 
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|A| 

= jo,  M = JO  + 1 

M — M — io  + 1 

€ 

a 

#Int 

m 

^25 

#Int 

m 

jmax 

#Int 

a 

3,4 

16,80 

8 

5,6 

108, 1264 

Q 

iwataiatM  % 

EE23H1 

a 

4 

18, 138 

8 

6 

142, 1603 

7 

224, 13168 

0.125 

a 

4,5 

26, 196 

El 

EH 

184, 2751 

o 

EH 

344, 18052 

0.0625 

a 

7 

281,3924 

8 

6 

442, 21924 

□ 

5,6 

44,454 

B9 

m 

396,8048 

11 

0.015625 

A. 

6 

i a 

8 

716, 12044 

m 

7 

i mm 

□ 

6,7 

82, 1318 

m 

8,9 

1152,27708 

m 

Tab.  1.  Estimated  and  determined  maximum  level  and  number  of  integrals  in  A*. 


kinds  of  effects  are  present,  namely  the  parametric  mapping  and  the  layer  near 
the  corner.  We  stress  that  we  do  not  intend  to  study  any  particular  choice 
for  the  adaptive  discretization,  i.e.,  the  choice  of  the  set  A.  We  are  primarily 
interested  in  the  behaviour  of  the  adaptive  approximation  c in  a realistic 
situation. 

As  trial  and  test  functions  we  used  the  biorthogonal  B-spline  wavelets 
on  the  interval  corresponding  to  the  parameters  d = d = 2 (i.e.,  piecewise 
linear  primal  functions  and  dual  functions  of  lowest  possible  order)  from  [11] 
(see  also  [7]  for  the  original  construction  on  1R).  For  the  system  0*,  we 
choose  as  in  [2]  piecewise  linear  interpolatory  wavelets.  This  of  course  implies 
that  the  computation  of  the  corresponding  wavelet  coefficients  cp  can  easily 
be  performed.  Moreover,  since  piecewise  linear  interpolatory  wavelets  are 
nothing  else  than  hierarchical  B-splines,  the  integrals  in  (10)  actually  only 
contain  scaling  functions. 

We  used  the  parameters  a = 2,  q = 3/4,  8k  ■=  (k  + l)-2  as  well  as  the 
corresponding  parameters  r — 2,  p = 3/4  and  8k  :=  k~L  for  the  right-hand 
side,  [2]. 

In  the  ID  tests  in  [2],  the  parameter  e was  chosen  as  the  error  in  the 
//'-norm  of  a corresponding  uniform  discretization.  From  a practical  point 
of  view,  this  is  of  course  unrealistic.  First  of  all,  the  solution  is  in  general 
not  known.  Moreover,  the  ultimate  goal  of  an  adaptive  scheme  is  to  avoid 
a (high  level)  uniform  discretization  but  to  use  the  degrees  of  freedom  in  a 
more  economical  way.  Hence,  we  performed  various  tests  on  the  choice  of  the 
parameter  e. 

Our  first  test  concerns  Je  in  Remark  2 and  the  number  of  integrals  needed 
for  computing  the  elements  of  the  stiffness  matrix.  Our  computations  are 
performed  in  this  way:  at  first,  we  start  from  the  minimum  level  (jo  = 3)  for 
the  used  wavelet  basis,  where  we  fix  a certain  e,  then  we  solve  the  problem 
with  scaling  functions  and  wavelets.  In  Table  1 we  compare  the  theoretical 
estimate  Je  on  the  maximum  level  in  A*  with  the  values  that  were  actually 
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determined  by  our  indicators  in  (13).  For  different  portions  of  the  stiffness 
matrix,  we  display  the  predicted  Je,  the  detected  maximum  levels  by  the 
indicator  as  well  as  the  minimum  and  maximum  number  of  integrals  needed 
for  computing  non  zero-entries. 

We  see  that  the  estimated  Je  is  always  larger  than  the  effectively  used 
maximum  level.  This  is  what  we  expected,  but  the  efficiency  of  the  algorithm 
may  be  reduced  by  an  excessive  over-estimate  of  Je . We  also  deduce  that  the 
efficiency  of  the  method  crucially  depends  on  the  choice  of  e since  the  number 
of  integrals  strongly  grows  for  decreasing  e.  This  is  surely  due  to  the  low  order 
of  the  interpolatory  wavelets  used. 

Next,  we  present  in  Table  2 the  average  number  of  integrals  computed 
for  the  stiffness  matrix  and  the  right-hand  side  for  the  first  two  levels  with 
the  same  e of  Table  1.  Here  J\  :=  max{|A|  : A £ A}. 


€ 

7.8e  - 3 

B 

ii 

CO 

4.04 

5.47 

hmi 

24.15 

-rt< 

II 

738.67 

/ 

9 

9 

HMEl 

II 

El 

Tab.  2.  Average  number  of  integrals  per  entry. 

We  deduce  that  the  choice  of  e not  only  influences  the  maximal  and 
minimal  number  of  integrals  as  shown  in  Table  1.  Since  the  average  number 
of  integrals  grows  when  e decreases,  the  choice  of  e effects  the  efficiency  of  the 
computation  of  the  whole  stiffness  matrix.  Moreover,  the  presence  of  the  first 
wavelet  level  also  increases  the  number  of  integrals. 

Finally,  we  consider  the  error  in  the  7/ 1 -norm  and  the  relative  error 


Mi,n 


for  different  choices  of  the  parameter  e.  In  Table  3 and  Figure  2 ‘rate’  cor- 
responds to  the  rate  of  convergence  w.r.t.  the  exact  solution  in  the  7/1  -norm 
for  the  first  two  levels  in  A.  We  see  that  these  quantities  do  not  depend  on  e. 
At  these  levels  the  relative  discretization  error  still  exceeds  the  relative  error 
(15).  This  explains  why  the  rate  of  convergence  is  basically  constant  w.r.t. 
the  choices  of  e. 

Hence,  the  choice  of  £ matters  only  if  this  value  is  at  least  of  the  same 
order  than  the  relative  discretization  error.  We  remark  that  also  for  increasing 
e all  scaling  functions  on  level  J whose  support  overlap  -t  belong  to  A* . 
This  implies  that  the  error  due  to  the  approximation  of  the  entries  of  the 
linear  system  is  bounded. 
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£ 

rate 

■Hi 

1.9666 

re?sEa 

11111 

0.0625 

0.03125 

gisrssi 

0.015625 

1.9591 

0.0078125 

1.9568 

0.1343 

Tab.  3.  Relative  error  and  rate  of  convergence  in  dependence  of  e. 


Fig.  2.  Rate  (left)  and  r^x  (right)  of  Table  3. 
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